

library(iterators)
library(foreach)
library(doParallel)
library(doMC)

library(plyr)

#################################################################
#### Set how many cores may be used for parallel processing. ####
registerDoMC(cores=detectCores())

###################################
#### Read in the matched data. ####

trad_gensurg <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_trad_patient_match_gensurg_11.2.2018.csv")
mod_gensurg <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_mod_patient_match_gensurg_11.2.2018.csv")
trad_ortho <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_trad_patient_match_ortho_11.8.2018_exact_procs.csv")
mod_ortho <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/simul_mod_patient_match_ortho_11.8.2018_exact_procs.csv")

trad_gensurg <- trad_gensurg[,which(!(names(trad_gensurg) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")))]
mod_gensurg <- mod_gensurg[,which(!(names(mod_gensurg) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")))]
trad_ortho <- trad_ortho[,which(!(names(trad_ortho) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")))]
mod_ortho <- mod_ortho[,which(!(names(mod_ortho) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")))]

##### Read in the readms_or_death_new variable. ####
tg_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_gensurg_new_readms.csv")
mg_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_gensurg_new_readms.csv")
to_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_ortho_new_readms.csv")
mo_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_ortho_new_readms.csv")

trad_gensurg <- merge(x=trad_gensurg, y=tg_readms, by=c("newID_1","newID_10"), all.x=T)
mod_gensurg <- merge(x=mod_gensurg, y=mg_readms, by=c("newID_1","newID_10"), all.x=T)
trad_ortho <- merge(x=trad_ortho, y=to_readms, by=c("newID_1","newID_10"), all.x=T)
mod_ortho <- merge(x=mod_ortho, y=mo_readms, by=c("newID_1","newID_10"), all.x=T)


###############################################
#### Read in the updated icu_yes variable. ####

trad_ortho_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/trad_ortho_new_cost_and_icu.csv")
mod_ortho_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/mod_ortho_new_cost_and_icu.csv")
trad_gensurg_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/trad_gensurg_new_cost_and_icu.csv")
mod_gensurg_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/mod_gensurg_new_cost_and_icu.csv")

trad_ortho_icu <- trad_ortho_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
mod_ortho_icu <- mod_ortho_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
trad_gensurg_icu <- trad_gensurg_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
mod_gensurg_icu <- mod_gensurg_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]


#####################################
#### Merge the new ICU variable. ####

trad_gensurg <- merge(x=trad_gensurg, y=trad_gensurg_icu, by=c("newID_1","newID_10"), all.x=T)
mod_gensurg <- merge(x=mod_gensurg, y=mod_gensurg_icu, by=c("newID_1","newID_10"), all.x=T)
trad_ortho <- merge(x=trad_ortho, y=trad_ortho_icu, by=c("newID_1","newID_10"), all.x=T)
mod_ortho <- merge(x=mod_ortho, y=mod_ortho_icu, by=c("newID_1","newID_10"), all.x=T)



gensurg_match <- rbind.fill(trad_gensurg, mod_gensurg)
ortho_match <- rbind.fill(trad_ortho, mod_ortho)

#### Modify the pair IDs to make them unique. ####
gensurg_match$pair_id_pats <- paste(gensurg_match$pair_id_pats,"g",gensurg_match$era_trad,sep="_")
ortho_match$pair_id_pats <- paste(ortho_match$pair_id_pats,"o",ortho_match$era_trad,sep="_")

#### Break each dataset apart by era. ####

trad_gensurg <- subset(gensurg_match, era_trad==1)
mod_gensurg <- subset(gensurg_match, era_mod==1)
trad_ortho <- subset(ortho_match, era_trad==1)
mod_ortho <- subset(ortho_match, era_mod==1)

both_match <- rbind.fill(gensurg_match, ortho_match)
trad_match <- subset(both_match, era_trad==1)
mod_match <- subset(both_match, era_mod==1)

#### Define the binary and continuous outcomes. ####

bin_outcomes <- c("readms_30","readms_or_death","dead_inhosp","ftr_30_ign_poa","ftr_30_use_poa","ftr_inhosp_ign_poa","ftr_inhosp_use_poa","icu_yes",
                  grep(x=names(gensurg_match),pattern="complic_inhosp_",value=TRUE), "dead_30","dead_60","dead_90","dead_180", "prolonged_d","step_up","complic_poa_any",
                 "readms30_post", "readms_or_death_new") # "robotic_assist_spx")

compsyn_pos <- match("complic_inhosp_COMPSYN", bin_outcomes)


bin_outcomes <- bin_outcomes[-compsyn_pos]

con_outcomes <- c("total_cost_complete","total_cost_index_only","anesth_time_claim","icu_days","los_long", "basic_los")

#################################################
#### Write functions to perform Gart's test. ####

garts_test <- function(data1, data2,variable, pair_id,treatment, comparison){
  
  #### Split the data apart by cases and controls. ####
  new_data1 <- data1[which(data1[,treatment]==1),c(pair_id,variable)]
  exp_data1 <- data1[which(data1[,treatment]==0),c(pair_id,variable)]
  new_data2 <- data2[which(data2[,treatment]==1),c(pair_id,variable)]
  exp_data2 <- data2[which(data2[,treatment]==0),c(pair_id,variable)]
  
  #### Rename the variables for convenience. ####
  colnames(new_data1)[2] <- c("var_t1")
  colnames(exp_data1)[2] <- c("var_t0")
  colnames(new_data2)[2] <- c("var_t1")
  colnames(exp_data2)[2] <- c("var_t0")
  
  #### Create new datasets that have 1 row for each patient pair. ####
  merge_data1 <- merge(x=new_data1, y=exp_data1, by=c(pair_id))
  merge_data2 <- merge(x=new_data2, y=exp_data2, by=c(pair_id))
  
  #### Get McNemar's table. ####
  mctable1 <- table(merge_data1$var_t0, merge_data1$var_t1)
  mctable2 <- table(merge_data2$var_t0, merge_data2$var_t1)
  
  #### Create a 2-by-2 table for Gart's test. ####
  gtable <- matrix(c(mctable1[2,1],mctable2[2,1],mctable1[1,2],mctable2[1,2]),nrow=2,ncol=2)
  
  #### Perform Gart's test, which is just Fisher's Exact test on the above matrix. ####
  gtest <- fisher.test(x=gtable, conf.int=TRUE, conf.level=0.95)
  
  
  #### Aggregate the results. ####
  results <- c(comparison,variable,gtest$method, gtest$estimate, gtest$conf.int[1], gtest$conf.int[2], gtest$p.value)
  names(results) <- c("Surgery", "Outcome", "Test", "OR Estimate", "OR 95% Lower Bound", "OR 95% Upper Bound","Fisher's Test p-value")
 
  return(results)
   
}

########################################################
#### Calculate Gart's test for each binary outcome. ####

registerDoMC(cores=detectCores())

garts_tests <- foreach(i=1:length(bin_outcomes),.combine='rbind') %dopar% {
  
    gt1 <- garts_test(data1=trad_gensurg, data2=mod_gensurg, variable=bin_outcomes[i], pair_id="pair_id_pats",treatment="new_surgeon_flag", comparison="gensurg")
    gt2 <- garts_test(data1=trad_ortho, data2=mod_ortho, variable=bin_outcomes[i], pair_id="pair_id_pats", treatment="new_surgeon_flag", comparison="ortho")
    gt3 <- garts_test(data1=trad_match, data2=mod_match, variable=bin_outcomes[i], pair_id="pair_id_pats",treatment="new_surgeon_flag",comparison="gensurg and ortho")
    gt <- rbind(gt1,gt2,gt3)
  
}

write.csv(x=garts_tests, file="/Users/sharpe/Desktop/Young Surgeons/table/garts_tests_REDONE_11.26.2018.csv")


#####################################################################################
#### Write a function to perform Wilcoxon Rank-Sum Test on continuous variables. ####


dnd_con <- function(data1, data2,variable, pair_id,treatment, comparison){
  
  #### Split the data apart by cases and controls. ####
  new_data1 <- data1[which(data1[,treatment]==1),c(pair_id,variable)]
  exp_data1 <- data1[which(data1[,treatment]==0),c(pair_id,variable)]
  new_data2 <- data2[which(data2[,treatment]==1),c(pair_id,variable)]
  exp_data2 <- data2[which(data2[,treatment]==0),c(pair_id,variable)]
  
  #### Rename the variables for convenience. ####
  colnames(new_data1)[2] <- c("var_t1")
  colnames(exp_data1)[2] <- c("var_t0")
  colnames(new_data2)[2] <- c("var_t1")
  colnames(exp_data2)[2] <- c("var_t0")
  
  #### Create new datasets that have 1 row for each patient pair. ####
  merge_data1 <- merge(x=new_data1, y=exp_data1, by=c(pair_id))
  merge_data2 <- merge(x=new_data2, y=exp_data2, by=c(pair_id))
  
  #### Calculate paired differences. ####
  merge_data1[,"difference"] <- merge_data1[,"var_t1"] - merge_data1[,"var_t0"] #### New minus Exp.
  merge_data2[,"difference"] <- merge_data2[,"var_t1"] - merge_data2[,"var_t0"]
  
  wt <- wilcox.test(x=merge_data2[,"difference"],y=merge_data1[,"difference"],correct=T, exact=F, conf.int=T,conf.level=0.95,mu=0)
  
  wt_vec <- c(comparison, variable, wt$method, wt$estimate, wt$conf.int[1],wt$conf.int[2], wt$p.value)
  names(wt_vec) <- c("Surgery","Outcome","Method","Difference Estimate","Lower 95% CI","Upper 95% CI", "P-Value")
  
  return(wt_vec)
  
}

continuous_vars <- foreach(j=1:length(con_outcomes),.combine='rbind') %dopar% {
  
  gensurg <- dnd_con(data1=trad_gensurg, data2=mod_gensurg, variable=con_outcomes[j], pair_id="pair_id_pats",treatment="new_surgeon_flag",comparison="gensurg")
  ortho <- dnd_con(data1=trad_ortho, data2=mod_ortho, variable=con_outcomes[j], pair_id="pair_id_pats",treatment="new_surgeon_flag",comparison="ortho")
  both <- dnd_con(data1=trad_match, data2=mod_match, variable=con_outcomes[j], pair_id="pair_id_pats",treatment="new_surgeon_flag",comparison="gensurg and ortho")
  all <- rbind(gensurg, ortho, both)
  
  
}

write.csv(x=continuous_vars, file="/Users/sharpe/Desktop/Young Surgeons/table/wilcoxon_rank_sum_tests_REDONE_11.26.2018.csv")





